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Abstract 

Conservation laws, in for example, electromagnetism, solid and fluid mechanics, allow an exact discrete representation 
in terms of line, surface and volume integrals. We develop high order interpolants, from any basis that is a partition of 
unity, that satisfy these integral relations exactly, at cell level. The resulting gradient, curl and divergence conforming 
spaces have the property that the conservation laws become completely independent of the basis functions. This means 
that the conservation laws are exactly satisfied even on curved meshes. As an example, we develop high order gradient, 
curl and divergence conforming spaces from NURBS - non uniform rational B-splines - and thereby generalize the 
compatible spaces of B-splines developed in We give several examples of 2D Stokes flow calculations which 
result, amongst others, in a point wise divergence free velocity field. 
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Be careful of the naive view that a physical law is a mathematical relation between previously defined 
quantities. The situation is, rather, that a certain mathematical structure represents a given physical 
structure. Burke ^ 

1. Introduction 

In deriving mathematical models for physical theories, we frequently start with analysis on finite dimensional 
geometric objects, like a control volume and its bounding surfaces. We assign global, 'measurable', quantities to 
these different geometric objects and set up balance statements. Take for example the global balance in ([1]), where 
the total mass/momentum/energy E inside a control volume V is only conserved (no change in time) if the in- and 
outgoing mass/momentum/energy fluxes Q over the bounding surfaces dV cancel. 

|£(y) = e(5y) — - |e = divq. 

(1) 

(discrete or global) (differential or local) 

This is exactly Gauss divergence theorem, depicted in Figure [1]:. Other balance equations in M? involve the funda- 
mental theorem of calculus, relating a global quantity associated with a curve L, to the values of a quantity at the 
boundary points dL, and Stokes circulation theorem, which relates the amount of rotation in a surface S to the amount 
of circulation around the bounding curve dS . 

While the association of physical quantities with geometry is clear in the global sense, it remains obscured when 
the mathematical model is written in local form, in ([TJ, as a differential equation. The local variables, i.e. the source 
field q and density e, obtained from a limiting process by shrinking the integration domain V up to a point P, although 
mathematically well defined, seem to have lost their geometric significance. 
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Figure 1 : The fundamental theorem of calculus. Stokes circulation theorem and Gauss divergence theorem. 



Classical numerical methods, in particular finite difference and nodal finite element methods, take the differential 
statement as a starting point for discretization. These methods expand their unknowns in terms of nodal interpola- 
tions only, and thereby disregard the underlying geometry of the physics. This can lead to instabilities, and perhaps 
more dangerously, to internal inconsistencies such as the violation of fundamental conservation principles. Where 
instabilities lead to outright failure of a numerical method, inconsistencies can lead to unphysical solutions that can 
go unnoticed by the human eye Js]]. This becomes more and more pronounced, as the trend is to simulate increasingly 
larger and complex non-linear phenomena, like multi-phase flows, fluid structure interaction and magnetohydrody- 
namics. 

To capture the behavior of a physical phenomena well, a discretization method should't only approximate the 
spaces of the infinite dimensional system, but should also follow the structure induced by the relations between them; 
in particular the structure induced by the fundamental balance equations depicted in Figure[T] By choosing degrees of 
freedom, not solely associated with nodes, but also with edges, faces and volumes in the mesh, we are able to exactly 
satisfy these relations in the discrete setting. 

We follow the pioneering work of Tonti flS, Mattiussi jill, Bossavit QHtl and the recent advances in Discrete 
Exterior Calculus [[ sji [Toll. F inite Element Exterior Calculus 111, 121. Compatible 111 [Tsl [l4i [Tsl [l6ll and Mimetic 
Methods \ VJ . 3 [l9[l2o[ 21 , 22 , 23 , 24 ] . These methods, do not focus on one particular physical problem, but identify 



and discretize the underlying structure that constitutes a wide variety of physical field theories. They are said to be 
'compatible' with the geometric structure of the underlying physics, i.e. they 'mimic' important properties of the 
physical system. This leads, amongst others, to naturally stable and consistent numerical schemes that have discrete 
conservation properties by construction and are applicable to a wide variety of physical theories. Furthermore, they 
offer insight into the properties of existing numerical schemes. 

We develop arbitrary order interpolants, starting from any basis that is a partition of unity, which satisfy the 
fundamental integral theorems exactly. The gradient, curl and divergence conforming spaces have the property that 
the conservation laws become completely independent of the basis functions. This means that the conservation laws 
are exactly satisfied at the coarsest level of discretization and on arbitrarily curved meshes. It is remarkable that 



inf-sup stability is automatically guaranteed when this physical structure is encoded in the discretization 112211231. 

As an example we apply our approach to NURBS (non-uniform rational B-splines) and thereby generalize the 
compatible spaces of B-splines introduced in yj. NURBS, the standard in Computer Aided Design (CAD), have only 
recently become popular in the finite element community due to the very succesfull IsoGeometric Analysis (IGA) 
paradigm ll25ll . IGA employs CAD technologies, such as B-splines and NURBS, directly within the finite element 
spaces, and thereby integrates computer aided design with finite element analysis (PEA). IsoGeometric Analysis not 
only bridges the gap between CAD and FEA, it also provides exact geometry representation at the coarsest level 
of discretization 12511 : possesses increased robustness and accuracy ll26ll : and refining strategies become practically 
applicable ll25ll . Development is ongoing and IsoGeometric Analysis has over the years reached a certain level of 
maturity. For an overview we refer the reader to ll27ll . 
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Figure 2: Physical quantities in M? can be associated witli eitlier inner or outer oriented points, curves, surfaces and volumes. They are horizontally 
connected by means of the balance equations, i.e. fundamental theorem of calculus. Stokes circulation theorem and Gauss divergence theorem, 
which ai'e metric free. The metric and material dependent relations, i.e. the constitutive equations represent the relation between physical quantities 
associated with dual geometric objects. 



1.1. Geometric structure of partial differential equations (PDE's) 

Geometry induces physics with a clear geometric structure. Consider Figure |2] Physical quantities are naturally 
related to geometric objects such as points, curves, surfaces and volumes. Furthermore, we can consider two separate 
types of orientation with respect to the geometric object: inner and outer orientation. Temperature, for example, is 
measured in a point; strain of a fiber along a curve; magnetic flux through a surface; mass flowing out of a volume; and 
an amount of rotation in a surface. Note there is room for interpretation since an amount of rotation could similarly 
be associated with a vortex filament. 

Physical quantities are horizontally connected by the fundamental theorem of calculus, Stokes circulation theorem 
and Gauss divergence theorem (Figure [TJ. These balance equations can exactly be described in terms of discrete 
physical quantities, that could have been obtained from some measurement process, and besides are independent of 
the shape of the geometry; i.e. they are intrinsically discrete and topological in nature. 

The metric, notions of length, area, volume, angle, along with material properties of the medium, are described 
by the constitutive equations, which require a differential formulation. These represent the relations between physical 
quantities associated with dual geometric objects. Consider for example Hooke's law which relates the strain of a fiber 
along an inner oriented curve, with the stress through an outer oriented surface. Another example is the proportional 
relationship between point-wise defined temperature and kinetic energy contained in a volume for a perfect gas. 

The schematic in Figure|2]can serve as a template to design compatible numerical methods for partial diff'erential 
equations (PDE's). By reformulating a PDF in terms of a set of first order equations one can distinguish between the 
topological nature of the balance laws and the metric nature of the constitutive equations. This idea is by no means 
new. In fact, Figure|2]is known as a Tonti diagram, after Enzo Tonti. For more details, we refer the reader to the work 
of Tonti |4l| and Mattiussi 10]. 

We consider Stokes flow in a domain Q, filled with an incompressible fluid with constant viscosity v. Under these 
assumptions. Stokes flow can be described by the following equations, 

vAu - grad p = f and div u = 0, (2) 

where u is the velocity, p the pressure and f the right hand side forcing. Using the the operator splitting, -Au - 
curl curl u + grad div u, and using the incompressibility constraint, div u = 0, we can decouple Stokes flow into the 
following first order equations, 

curl u - CO, curl (o + grad p = -f , div u = 0, (3) 
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Figure 3: Tonti diagram illustrating mixed formulation for Stokes flow with outer oriented variables. Velocity u is associated with outer oriented 
surfaces, pressure p is defined in outer oriented volumes and the vorticity <i) is related to outer oriented lines. 
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Figure 4: Tonti diagram illustrating mixed formulation for Stokes flow with inner oriented variables. Velocity u is associated with inner oriented 
lines, pressure p is associated with inner oriented points and the vorticity &) is related to inner oriented surfaces. 



where oi is the vorticity. Observe that Stokes circulation theorem (Figure[TJ)) means that u in curl u is associated with 
inner oriented curves, while Gauss divergence theorem (Figure [1]:) implies that u in div u is associated with outer 
oriented surfaces. The variable u thus has a different geometric interpretation depending on the type of balance law 
and needs to be treated accordingly. 

It is possible to discretize variables on staggered grids of opposite orientation, as in many finite volume methods, 
and make an explicit relation between quantities associated with dual geometric objects. We will however follow a 
more finite element type of method and circumvent the use of a dual (staggered) grid by using integration by parts. 
While grad, curl and div are the differential operators that naturally occur in the fundamental theorem of calculus. 
Stokes circulation theorem, and Gauss divergence theorem, we can define grad*, curl* and div* as their formal 
Hilbert adjoint up to a possible boundary term [21], 



(a, -grad* P) = (div a,p), (y, curl* 6) = (curl 7, 6), (e, -div* O = (grad e, O 



(4) 



While the grad, curl and div operators are purely topological (metric free) and allow an exact discrete representation 
(Section[3]l, the grad*, curl* and div* are metric dependent and can only be approximated in the discrete setting. 

We can use these operators in a mixed Galerkin setting, where the resulting mixed formulation depends solely 
on physical considerations. Take for example the mixed formulation of the Stokes problem depicted in Figure [3] 
where the variables are associated with outer oriented geometric elements. In this case we apply the topological 
equation div u = for conservation of mass, which will eventually give a point-wise divergence free velocity field 
in the discrete setting. Furthermore, inflow boundary conditions can be enforced strongly, while tangential velocity is 
prescribed in a weak sense. 

We could also choose to associate all variables with inner oriented geometric elements. This mixed formulation 
is shown in Figure |4] In this case the equation curl u = « is topological and allows an exact discrete representation. 
Furthermore, tangential velocity can be enforced strongly, while normal velocity is prescribed in a weak sense. 

The paper is structured as follows. In Section |2] we introduce differential forms as the natural mathematical 
objects representing physical fields. Differential forms allow a continuous description of physics that maintains its 
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geometric content and allows us to separate the topological from the metric dependent structure. Topological structure 
is intrinsically discrete, and can exactly be encoded in terms of discrete objects known as chains (discrete geometric 
objects) and cochains (discrete quantities), see Section |3] Metric structure, on the other hand, requires a continuous 
description. In Section |4] we derive a new procedure to develop high order cochain interpolants, from any basis that 
is a partition of unity, that by construction satisfy the balance laws. As an example, we develop high order gradient, 
curl and divergence conforming basis functions from NURBS and apply these newly developed spaces in a mixed 
Galerkin setting, Section|5j to Stokes flow and perform some numerical calculations to asses the numerical approach 
(section|6]l; finally, conclusions are drawn in section]?] 



2. Differential modeling 

The mathematical description of physical phenomena such as electromagnetism, solid and fluid mechanics, relies 
heavily on the use of line, surface and volume integrals. The differential objects that appear as their integrands 
are called differential forms (see Example 12.1b . and are studied in the mathematical field of differential geometry 



A description of physics in terms of differential forms offers significant benefits. Most notable, differential forms 
maintain a clear relation with the underlying geometry and therefore allow a separation of the metric dependent 
content, of a physical theory, from its topological (metric-free) part. This relation with geometry also makes it possible 
to transfer operations on the geometry to operations on the variables associated with that geometry. Furthermore, 
representation of variables in terms of differential forms offers a generalized concept of derivative, and thereby reduces 
the fundamental theorem of calculus. Stokes circulation theorem and Gauss divergence theorem to one equation, 
known as the generalized Stokes theorem. Most material presented in this section can be found in the references cited 
above. 

2.1. Differential forms 

Consider a sufficiently smooth «-dimensional domain D. with boundary dQ. and local coordinates x = (x' , x^, x"^ 
A A:-form a*" is a mathematical expression of the following form, 

a'' ai(x) dx' with dx' - dx'' A dx'- A ... A djc", (5) 

where the indices satisfy, 1 < ii < (2 < ... < it < n, the fl/(x) are smooth functions that denote the spatial distribution 
of some physical quantity, and the basis elements dx' refer to the associated geometry. The collection of all fc-forms 
is a vector space A*(Q) of dimension equal to the binomial coefficient Some examples of differential forms in 
are depicted below. 

Example 2.1 (Differential forms in M?). Temperature, T^, is a 0-form, a scalar function that assigns to every point 
X in domain Q. a real valued temperature; force, F^, is a 1-form, since it can naturally be integrated along a curve to 
obtain work; flux density, Q^, is a 2-form and can be integrated over a surface to obtain the global flux; mass density 
p^, is a 3-form, since it is readily integrated over a volume to yield mass. 

0- form Temperature: T'^ - T(x), 

1- form Force density F^ - fi{x)dx^ + f2{x)dx^ + f2(x)dx^ , 

2- form Flux density Q - qi(x)dx A dx' + q2{x)dx A dx + ^3(x)dA A dx , 

3- form Mass density - p(x)dx' A dx^ A dx^*. 



The A, is called the wedge product, which is a map, A : A*(f2) x A'(f2) i-> A''^'{Q.), k + l < n. The wedge product 
is skew symmetric, a'' A b' = (-1)*^^' b' A a''. So, dx' A dx^ - -dx-' A dx', which implies that the orientation changes 
under a permutation of the elements and Unear independence, since dx' A dx' - -dx' A dx' - 0. 
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A;-forms are the natural integrands over A:-dimensional geometric objects, 

(fl^Qi):= r fl*€M (6) 

JSit 

Here (., .) denotes duality pairing between the ^-form and a ^-dimensional geometric object. This duality, between 
geometry and physics has important consequences. It induces physics with a clear geometric structure. This allows us 
to separate the metric and material dependent content in a physical theory from the topological, metric independent 
part. Furthermore, it makes it possible to transfer operations on the geometry to operations on the variables associated 
with that geometry. 



2.2. Topological structure 

Integration itself is a metric free operation. If €) : f2' i-> Q, is a smooth map from an n-dimensional reference 
domain Q' = [0, 1]" to the «-dimensional physical domain Q, the pullback, O* : A*(f2) i-> A'^(Q.'), maps A:-forms on 
Q to A;-forms on Q.'. This means we can perform the integration of a fc-form in the reference domain Q', 

r r ^ (fl^(l)(Q^)) = ((DV,Q[). (7) 

J*(n;.) Jn;. 

This means that the pullback <t)* is the formal adjoint of the map O in the duality pairing defined in (|6]). Important 
properties of the pull back are linearity and commutation with the wedge product, O* (a'' A b'^ - (^*a''^ A ^O*^'). 

Topological structure is induced by balance laws, which relate a quantity associated with a geometric object to 
another quantity which is associated with its boundary. In the theory of differential forms, this balance is represented 
by the generalized Stokes theorem, the mother of all equations, 

r dw*-i = r w*^-' ^ (dw*-',ni.) = (8) 

Again we can observe the duality between geometry and physics: the boundary operator, d, is the adjoint of the 
exterior derivative d. While the boundary operator is a map d : ^k-i, the exterior derivative is a map, 

d : A^-i(f^) ^ A^(f^)- 

The exterior derivative d is a coordinate independent generahzation of the well known vector calculus identities, 
the gradient, curl and divergence operators in R"*. The Stokes theorem ([8]) thus generalizes the fundamental theorem 
of calculus {k = 1), Stokes circulation theorem (k - 2) and Gauss divergence theorem [k = 3), depicted in Figure 
[1] Furthermore, the exterior derivative is no more difficult to compute in a curved coordinate system than it is in a 
Cartesian frame. One simply takes the differential of the components and follows the properties of the wedge product. 



Example 2.2 (Action of the exterior derivative in M^). Consider local coordinates (x' , x-, jc^). The action of the 
exterior derivative on the 0-form T^, 1-form and 2-form from ExamDle \2.1\ is given by 

. n dT , dT , J dT . 
dr° = — dx' + —dx" + —dx' 

ox^ ax^ ox^ 

Observe that the components are given by those of the familiar grad, curl and div operators, respectively. 
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The exterior derivative has a number of important properties: it is a Hnear operator (|9]i; satisfies a Leibniz rule for 
differentiation (fTOl i: its metric free nature is reflected by commutation with the pull back (fTTT l: and is exact ( fTSl i: 

d (fl^ + b'') = dfl* + db'^ (9) 

d(fl'^ A b') = dfl*^ A - (-l/fl* A db' (10) 

d(D* = (t*d (11) 

dod = (12) 

This last property is analogous to the vector calculus identities, curl grad - and div curl - 0. 

The n + 1 -spaces of differential forms in an «-dimensional domain Q., satisfy the following sequence, known as 
the de Rahm complex, 

R A^Q) — ^ A'(Q) — ► ... ^ A"(Q) — " 0. (13) 



which is exact on contractible domains. The main focus of this paper, see SectionH] is to construct discrete spaces of 
differential forms, that follow the same structure, such that conservation and balance laws can be strongly enforced. 



2.3. Metric structure 

Constitutive equations depend on the local metric, notions of length, angle, area, volume etc., and material prop- 
erties of the medium under consideration. In order to measure the local metric we require a point wise inner product 
of A;-forms, (., .) : A* (Q) x A*^ (Q) i-> R. In particular, a Riemanian metric gives rise to the Hilbert space inner 
product on A'^(Q), 

(a\/?*)^:= r (Q'*,/?*)dO= f a^A-kp^. (14) 
Here ★ denotes the Hodge star operator, which is a map ★ : A*^ (Q) i-> A" (f2). 

Example 2.3 (Action of the Hodge star in M^). In with orthononnal coordinates y - {y^,y^,y^^ we have 
* 1 = dy' A dy~ A dy^ , ★ dy' A dy~ A dy"' = 1 



★ dy' = d/ A dy^ -k dy~ - dy^ h dy^ , -k dy^ - dy^ A dy^ . 

The Hodge not only maps A:-forms into (« - A:)-forms, but also changes its type of orientation from inner to outer 
and vice versa. We can therefore apply it to connect two copies of an exact sequence into the following structure, 

M A°(Q) — ^ K\Q) — ^ — ^ A^(Q) 

* * * * (15) 

a3(0) ^ — A2(Q) ^ — A'(Q) ^ — A°(n) ^ M 

Note that this is exactly the structure given in Figure|2l While the exterior derivative d models the topological structure 
of the balance laws, the Hodge ★ models the metric structure of the constitutive equations. 

Using Leibniz's rule for differentiation ( fTOt . and noting that a double application of the Hodge star changes nothing 
up to a possible minus sign, ★★ = (-1)'^'^""'^^ we can derive an adjoint operator for the exterior derivative, 

d(a*-' A = do-*-' A *y6* - (-l)*^*-' Ad*/?* 

= da*-' A ★/?* - (-l)"(*+i>+io'*-i A * (★d*)/?*. 
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Figure 5: The discrete Hodge invokes a connection between variables associated with dual geometric objects. 



where d* = (_i)n(*:+i)+i is the codifFerential. Using the inner product property of the Hodge (fT4l) and using 

integration by parts, we find that the codifFerential is the Hilbert adjoint of the exterior derivative with an additional 
boundary term, 

(da^-',/3*) - (a'^-',dV) = f a^-' A (16) 

The coderivative d* is thus a generalization of the grad*, curl* and div* introduced in (jUi. Equation ( fT6l l will prove 
to be very important in the mixed Galerkin setting explained in Section |5] The exterior derivative and codifFerential 
can be combined to construct the Laplacian A. The Hodge Laplacian is a map, A : A* (Q) A*^ (Q), 

A = dd*+d*d. (17) 



3. Discrete modeling 

Classical numerical methods, in particular finite difference and nodal finite element methods, expand their un- 
knowns in terms of nodal interpolations and run into trouble when it comes to conservation. Conservation, by defini- 
tion (Generalized Stokes theorem (|8]l) is a relation between a global 'measurable' quantity associated with a geometric 
object and another global 'measurable' quantity associated with its boundary. By choosing degrees of freedom asso- 
ciated with nodes, as well as edges, faces and volumes in the mesh, we are able to exactly satisfy these relations in the 
discrete case. 

The constitutive equations describe a relation between physical quantities associated with dual geometric objects. 
In discrete space, this is modeled by a discrete Hodge * operator which invokes a relation between staggered grids of 
opposite orientation, see Figure |5] Staggered finite volume methods explicitly build such a dual (staggered) grid and 
thus explicitly construct a discrete Hodge to connect variables on both grids. In this paper we use a Galerkin finite 
element type of approach and circumvent the use of a dual grid by applyin g integr ation by parts, using ( fTSI l. 



In this section we introduce geometry in terms of algebraic topology OOl bill . Topology describes the relations 
between oriented geometric objects (chains) and what 'lives' on those objects (cochains), however, without the notion 
of distance or measure. Algebraic topology thus allows us to encode the purely discrete and topological content of the 
physics into the discrete model, without any approximation. 

Algebraic topology can in many ways be regarded as the discrete analogous of differential geometry. The duality 
pairing between chains and cochains is analogue to that of integration of differential forms. By defining a formal 
adjoint of the boundary of a chain, we obtain a discrete derivative acting on cochains, known as the coboundary 
operator. This discrete derivative is constructed such that it exactly satisfies Stokes theorem in terms of chains and 
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(a) Invertible map (b) Numbering and orientation 




(c) Partitioning of reference domain in terms of of 0-cells (points), 1-ceIls (edges), 2-ceIls (faces) and 3-cells (volumes) 

Figure 6: (a) The curved domain <I)(n') = n has the same topology as the reference domain n'. (c) Once the reference domain Q' is partitioned 
into 0-cells (points), 1-cells (edges), 2-cells (faces) and 3-cells (volumes), every cell can be given an arbitrary numbering and orientation (b). 

cochains. Because all objects and operators that will be introduced in this section are topological and thus completely 
metric-free, they have the same form and value on topologically equivalent grids. The discretization of these relations 
is equivalent on a nice uniform Cartesian mesh as it is on a highly curved grid. Furthermore, they do not change on 
moving meshes as long as the topology stays the same. 



3.1. Cell complexes, chains and the boundary operator 

The pairing between physical quantities and oriented geometric objects, such as points, edges, faces and volumes, 
leads to a straightforward discretization of the «-dimensional computational domain Q. The topological description 
of Q. is given by a so called cell complex D, a formalized concept of the discretization of space, which is a partitioning 
in terms of A:-dimensional sub-domains {k - 0, n) called A;-cells, which we denote by cr;^) 

Figure |6] depicts a cell complex in M?, where we can distinguish between 0-cells (points), 1-cells (edges), 2-cells 
(faces) and 3-cells (volumes). For the present work it is most convenient to think of a cell complex as a union of 



the topological viewpoint, both descriptions are equivalent 1132 1. 



implex as 

^-cubes. An alternative would be to partition the cell complex in A:-simplices as is done in e.g. IglllOl ll,[12l]. From 



Once all A:-cells in the cell complex D have been numbered and given an orientation (see Figure [SJ)), they can be 
collected to form a ^-chain, aj^) e C{k){D), 



C(k) = [c 



-(^(kulj^,, withc^e {-1,0,1). (18) 



In the representation of a chain we use superscript for the coefficients c-' and subscript for the basis A:-cells crj^j y. 
Chains can be used as a discrete representation of the geometry. In the description of geometry, however, we will only 



be concerned with chains with coefficients, d of {-1,0, 1). These correspond to either a cell with orientation opposite 
to the chosen default one, a cell not part of the chain, or a cell with default orientation. In actual matrix calculations 

we directly use the coefficients as the column vector c -{cP c' • • ■ c^) . 

Connectivity between volumes, faces, edges and points, is encoded in the boundary operator, d. The boundary 
operator is a linear map d ; Qa) i-> C{k-\) defined as. 



dci 



ik) 



The boundary 5 of a A:-chain returns a unique {k — l)-chain and can be calculated as a linear combination of its A:-cells 
boundaries. More precisely, do-(^k).i - [e'j cr(A.-i),,j.^, where e^ has the value, 

1. e'j - 0, if cr(A_i) , is not part of the boundary of cr(A) j. 

2. e'j - 1, if cr(A_i) , is part of the boundary of cr^ii-, j with compatible orientation. 

3. = -1, if cr(A_i)_/ is part of the boundary of o-(k)j with incompatible orientation 
Figure|2]illustrated what we mean with compatible and incompatible orientation. 






Figure 7: Compatible orientation (+) and incompatible orientation (-) between volumes and boundary faces, faces and boundary edges, edges and 
boundary points. 

Hence, 



Xk) 



So the boundary of C(A) is the unique (k - l)-chain, c/(a-i), of which the coefficients are given by d' = Zij=o'"'^}- 
The action of the boundary operator therefore allows the matrix vector product d = Ei_i,iC. Here Ei_i i is a 
rank(C(<:_i)) x rank(C(A.)) incidence matrix with coefficients (E/;_i jt),-^ = e'j. Note that in matrix calculations we only 
use the coefficients, not the basis chains. Example l3. ll shows the incidence matrices associated with the numbered and 
oriented cell complex of Figure |6j). 

Example 3.1 (Incidence matrices in M.^}. Consider the numbered and oriented cell complex in Figure^jp. The con- 
nectivity between points, edges, faces and volumes is encoded in the following incidence matrices. Eq.i maps from 
1-chains to 0-chains, Ei _2 maps from 2-chains to 1-chains ant/E2,3 maps from 3-chains to 2-chains. 
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Figure 8: A double application of the boundary operator leads to an empty chain. Note that all edges are oriented in the opposite way, so any value 
associated with them cancels. 

One can readily check that the boundary of the boundary is empty, 

Eo,iEi,2 = (O ... O)^ and Ei,2E2,3 = (O ... O)^ 

We now introduce the most important property of the boundary operator, which leads to some profound conse- 
quences in the description of geometry and physics. Picture for example a finite dimensional volume V in M^. It is 
clear that its bounding surface dV encloses all of V (this is in fact the definition of the boundary). This means that the 
surface dV has no boundary itself; it is boundaryless. Figure[8]illustrates that taking the boundary twice of a volume 
leads to an empty 1 -chain. 

In general, taking the boundary twice of a ^-chain leads to an empty {k - 2)-chain, 

ddc(k) = 0(^2) for all q^., € C(^)(D). (19) 

The set of A:-chains and boundary operators thus gives rise to an exact sequence on contractible domains, the chain 
complex (C(j:)(D), d) 

■ ■ ■ " Ca-im ^ Qk)(D) ^— C(,+i)(Z)) ^ (20) 

3.2. Cochains and the coboundary operator 

So far we have learned how to calculate with discrete geometric objects, such as oriented points, edges, faces 
and volumes. By duality pairing with discrete geometry, we introduce cochains as discrete analogues of differential 
forms. Cochains can describe for instance integral values along a chain. Integration itself is however nothing but the 
'measurement technique' that assigns global values to cells. For our purpose, it is most convenient to think of cochains 
in the most general way as any set of global values associated with oriented points, edges, faces and volumes in the 
mesh (cell complex). More clearly, they need not be integral values. This fact will allow us maximum freedom in the 
projection of differential forms - to be discussed in the next section - where we use cochains as degrees of freedom in 
combination with suitable basis functions based on NURBS. 

A A:-cochain is an expression which looks like, 

fl*** = {fl; ■ o-'*>''}'^y withfl, eM and o-**>''' € C**>. (21) 

In contrast to chains, we use subscript for the coefficients a, and superscript to denote the basis cr^*^' '. While we 
considered chains only with coefficients {-1,0, 1), we allow the coefficients of cochains to be any real number As 

with chains, we shall use the coefficients directly as the column vector a = (flo ai ■ ■ fl ?) in matrix calculations. 

Cochains are linear functionals on chains, and by choosing the basis cochains cr**^*"' dual to that of the space of 
chains, l^cr^''^-' , cr^/i) Jj = we can define the duality pairing between chains and cochains as, 

■S' 

C(t)) := ^ fl/ • c' = a^c e E (22) 

1=0 
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Figure 9: A double application of the coboundary operator leads to an empty cochain. Note that the two additions of face values per volume cancel 
each other out. 



Note the similarity between duality pairing of differential forms and geometry by means of integration In the 
discrete setting, integration is replaced by summation (l22l i. 

We can now define a discrete Stokes theorem, the mother of all equations, in terms of chains and cochains. As in 
the continuous setting where the exterior derivative is the formal adjoint of the boundary operator ([8]l, we can define 
the coboundary operator on cochains as the formal adjoint of the boundary operator on chains, 

(23) 

for all i)**^"'^ 6 C(k-\){D) and C(t) € C(k){D). While the boundary is a map, d : Qjt) i-> C«_i), the coboundary is a map, 
6 : C**"'' i-> C**^'. Analogous to the exterior derivative d acting on forms in the continuous setting, the coboundary 
operator acts as a discrete derivative on cochains. The coboundary is thus a discrete version of the gradient {k = 1), 
curl (k - 2) and divergence {k - 3) operators from vector calculus. 

Earlier we remarked the fact that d o d - , has important consequences in physics. Amongst others, it implies 
that 6 o 6 = 0, see Figure |9] This is easily proven using ( l23T l twice and fT% . 

This property is analogous to (fT2l i and the important vector calculus identities curl grad = and div curl = in 
continuous space. 

We can now set up the following exact sequence, known as a cochain complex (c^''\D), S^, 

C(K-i\D) C^k)^D) C^'^'HD) . ■ ■ • (24) 



Since we have a matrix representation of the boundary operator, we can derive a matrix representation of the 
coboundary operator, using the adjoint property of the boundary and coboundary operator ( l23l l. 

(M*-", 5c(,)) P (E,_u.c) = (E[_,,,b)' c ^ {6b^'-'\ c(,)) . (25) 

We can conclude that the matrix representation of the coboundary operator is given by D^ ^-i = Ej^ ^. We can 
therefore perform difi'erentiation in the discrete setting using the matrix representation a - D^ ^-ib. Furthermore, 

^k+i,kDk.k-i - (p ■ ■ ■ O) , SO the matrix representation of the gradient, curl and divergence exactly preserves the 
null space of the differential operators. 



4. Commuting projection 

Balance equations, represented by the Generalized Stokes theorem ([8]l, allow an exact discrete representation in 
terms of chains and cochains using ( l23b . Since the structure of these equations is inherently discrete and metric free, 
we do not expect that basis functions play any role here. The constitutive equations - the material and metric dependent 
relations - on the other hand, require a continuous formulation of both the geometry and the field variables, and here 
is where the basis functions come in to play. 
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A continuous representation of field variables involves approximation by a suitable projection, tt/, : A*(Q) i-> 
A^j(f2), from an infinite dimensional space A'^(n) of differential forms, to a finite dimensional conforming subspace 
A^j(f2) c A'^(Q). Approximation always involves a loss in information. It is however important that the error is 
bounded by a certain constant C < oo. 



approximation property: 



for all fl* € A'-(Q) 



(26) 



for a norm defined on A*(Q), where p denotes the approximation order and h is a measure of the the maximum 
partitioning size. In order to be consistent, the projection should reproduce all of A*^(f2), which means that, 

for all fl^ e A*(Q) 



consistency property: 



k k 

TTl, O TT/, fl — TTIjO 



(27) 



These two requirements are in general met by polynomial approximations. For interpolation estimates of for example 
NURBS and splines, see l33ll34ll respectively. 

To capture the behavior of a physical phenomena well, a discretization method should not only approximate 
the spaces of the infinite dimensional system, but should also follow the structure induced by the relations between 
them, i.e. the structure induced by the balance laws and constitutive equations (Figure |2]i. Only then will a discrete 
representation of the physics have any physical significance. The focus of this paper is to derive a projection of 
differential forms which commutes with diff'erentiation. 



A* (Q) ■ 



K (") ■ 



A*+i (Q) 



(28) 



(Q) 



Projection and subsequent differentiation should give the same result as first applying the derivative and then do 
the projection. This commuting diagram property will guarantee that conservation and balance laws remain exactly 
satisfied in the discrete setting. Furthermore, this property is the key to naturally stable and consistent numerical 
algorithms. It can be shown that ( l28l l implies a conforming Hodge decomposition, which together with the Poincare 
inequality proves inf-sup stability. For more details we refer the reader to the papers by Kreeft et. al. ll22ll23[l . 

To assure that the balance laws are exactly satisfied, even on curved meshes, we require that the projection is 
independent of geometric transformations. This means that the projection should be constructed such that it commutes 
with the pull back as well. 



A* (<1)(Q')) 



Aj,(0(Q')) 



■ A*^ (no 



(29) 



■Af,(Q') 



This can be achieved by constructing the projection such that it preserves the integrals of a differential A:-form, a*^, 
along a chain C(ic). Using the fact that the pull back (D* is the adjoint of the map <1> we obtain the desired result. 



(30) 



4.1. Reduction, change of bases and reconstruction 

We extend the commuting diagram in (l28T l to define operations between the continuous (Section |2]i and discrete 
formalism (Section [3]l. Following Bochev and Hyman lIlSll . we introduce two separate operators, the reduction K 
and reconstruction I - which map differential forms to cochains and cochains to a finite dimensional representation 
of differential forms, respectively. Additionally we introduce a change of basis from cochains to a different type of 
cochains. As we shall see, this change of basis will allow maximum freedom in the choice of reconstruction method. 
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Reduction "R : A*(n) C<*'(D) is an abstraction of the measurement process and is given by the DeRahm map 
(|3TT) . which is defined by means of integration, 



(31) 



To be able to use the coboundary as a discrete derivative, as outlined in Section [3] we require that the reduction 
commutes with differentiation. 



A*(n) ^ A*^+i(Q) 



C(k) (D) — ^ C<*+'> (D) 



(32) 



which is indeed the case as can be proven by using Stokes theorem and the duahty between the boundary and cobound- 
ary operator, 

{K {da') , q,,,,) ^ f da* i r a'^{<R («*) , dc,,,,,) P {Sli {a') , q,,,,) 



Change of basis T : C'-''\D) i-> C*^'\D) is an invertible map from a cochain a**' - {a^ ■ cr**''-'} '_„ ^ C^'^\D), to a 



different type of cochain a**' = \aj 



cr 



)/■=() 



e C^''\D) which are defined by means of integration ( ISTT i. The change 
of basis makes it possible to use a much wider range of reconstruction methods. It provides the freedom to use 
degrees of freedom other than nodal values, curve, surface and volume integrals, while maintaining all the advantages 
of a compatible approach. For the reconstruction of 0-forms, for example, we no longer require nodal interpolants, 
but can use any basis that possesses partition of unity, for example NURBS. In practice, the change of basis is given 
by a square invertible matrix equivalent to that in an interpolation / histopolation problem. 

Since 'F is a map from cochains to cochains, we require that it commutes with discrete differentiation by means 
of the coboundary operator. 



C^*-') (D) ■ 
r- 
C^*-') (D) ■ 



■ C(*^+i> (D) 

■ C(*^+'> (D) 



(33) 



The commuting diagram in ( l33l l guarantees that the mathematical structure remains exactly the same after change of 
basis. The discrete modeling tools introduced in Section [3] therefore remain unaltered and we can perform discrete 
differentiation using the matrix representation Dj^+ijt = Ef, . , as explained in the end of Section[3] 



Reconstruction I : C^''\D) ^ A*(Q) maps cochains a^'^' - [aj ■ c^*'"'' } .^^ ^ C^'^\D) back to a finite dimensional 
representation of a differential form ajj e Aj^(Q) c A'^(Q). In order to use the discrete calculus introduced in Section 



[3]the reconstruction is required to commute with differentiation. 



C(*) (D) ■ 



(Q) ■ 



■ C('^+'> (D) 



(34) 



This way we can perform the gradient, curl and divergence in the discrete setting, while continuous representations 
can be reconstructed where and whenever required. 
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The commuting property of projection with differentiation (l28T l can be seen as the composition of the commuting 
relations in ( l32b . ( 1331 ) and ( l34b and is illustrated in the following diagram. 



A^' (Q) A*+i (Q) 



(35) 



The choice of reconstruction, I : C'-'^HD) A* (Q), uniquely defines that of J : C<*+"(D) A*+' (Q) 
through commutation with differentiation ( l34l i. This in practice means that we choose a basis for A'| (Q) and derive 
the remaining finite dimensional spaces of differential forms, AjJ (Q), for ^ = 1, «, using the commuting properties. 

Furthermore, from ( IZTT i it naturally follows that the change of basis is related to the reduction and reconstruction 
by, 'F := Hoi. Therefore the choice of basis for 0-forms determines not only the basis for all other spaces of 
differential forms, but also uniquely defines the change of basis !F : C^*'(D) i-> C***(D) for A: = 0, 1, «. 

We will study the commuting properties in ( l35l l in more detail in the univariate case, see Figure [10] Once all 
spaces and operators are determined in the univariate case, we shall use the tensor product to develop multivariate 
spaces of discrete differential forms. 



4.2. Interpolation and histopolation 

The commuting relations in (l35l l are graphically illustrated in Figure [TOlin the case of ID-space. Projection of a 
0-form is equivalent to solving an interpolation problem where a set of data measurements are interpolated, while the 
projection of a 1-form is equivalent to solving a histopolation problem, where a set of line integrals are preserved in 
the projection. 

Say we wish to approximate the 0-form in Figure[TOl which represents for example a temperature field T{x). We 
seek the finite dimensional approximation Th{x) - n^Tix) in the following space, 

a0(Q) |j [Tj ■ ^r*°'-^}"^^, = J Tj Njix), for all Tj 6 k| , (36) 

where the degrees of freedom Tj are the coefficients in the 0-cochain jfy ■ a"*'^''-'| _^ and the |A^y(jc)| ,^ are linear 

independent basis functions that posses a partition of unity, Yj'j=o ^/-''^) - 1 ■ 

The projection of a 0-form Ti,(x) - 7Ti,T{x) is required to follow K {T{x)) - H {Ti,{x)), which means that Ti,{x) 
interpolates T(x) at a chosen set of nodes {-f,)"^,). Then, using (|36] |. we can set up the n + 1 by n + 1 system of linear 
equations, 

n 

nhT{xd = Yj Njixi) = T(xi) for / = 0, 1 , n (37) 

j=o 

and find the unique set of coefficients jfyj,^. This system is guaranteed to have a solution when the matrix N with 
coefficients N,^ = Njixj) is invertible, which is the case when x,- e span {A^,(x)), for / = 0,...,n Issll . Note that the 
change of basis T -Ho I thus given by this square interpolation matrix N. The inverse map ;F ' - N depicted 
in the middle of the left hand side of Figure [TO] is thus equivalent to solving the interpolation problem. 
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Figure 10: While a zero form is projected by interpolating a set of nodal values (left), its deiivative (light) is projected by interpolating a set of 
integral values; a process called histopolation. This construction guarantees that the projection commutes with diiferentiation. Furthermore, the 
projection is broken down in three stages: 1) the reduction process K reducing a 0-form and its derivative to a set of nodal and integral values, 
respectively; 2) A change of basis from nodal and integral values to new 'node' and 'edge' type of degrees of freedom. This change of basis is 
equivalent to solving an interpolation problem N"' and histopolation problem M"', respectively. 3) the reconstruction process I which is simply 
given by a linear combination of the new 'node' (0-cochain), or 'edge' ( 1-cochain) degrees of freedom with the appropriate basis functions [Ni(x)]"_Q 
and {M,(x)|"_j respectively. The resolution is purposely kept low to illustrate the concepts of reduction, change of basis and reconstruction. 
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If we would like to approximate the derivative of the temperature field, the 1-form m'(x) = dT{x), we require a 
basis for 1 -forms, 



A,^(Q) = \l[uj ■ <5-'"'^]|'^j = ^ Uj Mj{x), for all e M i (38) 



where the degrees of freedom uj are the coefficients in the 1-cochain ^uj ■ cr*'''-' j _^ and therefore the basis functions 

|My(x)| can be associated with edges. 

The projection of a 1-form m'^ = ;t/,m' is required to follow ^m') - H {u^^^, which means that the approximation 
should preserve the line integrals of the 1-form u^{x) - d T(x) along the edges e, = [x,_i, jc,]. We can therefore set up 
the following « by n system of linear equations, 

r ttam'^^m; r M/a)= r M'(x) = r(x,)-r(x,_i) for/=l,2, (39) 

This process (right side of Figure [TOl i is called histopolation, from interpolation of a histogram. In this case, the 
change of basis 'F -"R o I,is equal to the invertible histopolation matrix M, of which the coefficients are given by, 

Mij^ j^ Mj(x) 

The inverse map ;F ' - M depicted in the middle of the right hand side of Figure [10] is thus equivalent to solving 
the histopolation problem. 

4.3. Edge functions 

The idea now is to choose a basis {A^,(x))"^q for 0-forms, and to derive the new edge type of basis functions 
{Mi{x)]"^i by using the commuting properties in (|34] |. Connectivity of edges e, and points x,-, c?e, = x, - x,_i, implies 
that we can apply discrete differentiation using the coboundary ( l23l l, m, - Ti - f,-! , for i - 1,2, ...,n. Hence, every 
coefficient f ,■ can be written as f,- - To + Yj'j^i {Tj - Tj-ifj = Tq + Yj'j^i 

Using this result in the projection of a zero form and applying the exterior derivative. 



dm, T{x) = d ^ f Ni(x) = f d ^ Ni(x) + ^ j^uj 



dNiix). 

/=0 \j=l } 



Partition of unity implies that, d 2"=o ^ - 0; hence the first term cancels. The second term can be rewritten such 
that we obtain a relation for the edge basis functions {M,(x))"^j, 

duh r°(x) = Ml dNx{x) + 

n n 

= {dA^i(x) + ... + dNnix)] Ml -H \dNi(x) + ... + dA?„(x)) M2 + - = ^ ^ dNj{x). 

i=\ j=i 

By defining the edge basis functions {M,(x))"^j as, 

n /—I 

Mi{x) = dNjix) = - ^ dNjix), for / = 1 , 2, n, (40) 
we have obtained the commuting projection, 

n n n 

1=0 !=1 1=1 
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(a) 



(b) 




(c) (d) 

Figure 11: Cubic NURBS basis functions (left) and associated quadratic edge functions (right) corresponding to a knot vector t = 
{ 0,0,0,0, 1,2,3,4,4,4,4 } under the influence of a dift'erent choice in weights. Every weight w, associated with basis function Ni(x) is set 
equal to one, except W3: Figure (a) and (b) correspond to 1^3 = 1 and thus simplify to B-splines; (c) and (d) correspond to 1^3 = 1/2. Note that the 
'node' NURBS basis features partition of unity 2"=o ^ii-^) = 1 on [0,4], while the edge functions feature unit integral Mj(x) = 1. 



As constructed, we can perform differentiation entirely independent from the basis functions. In fact, in dr" — 
Ij''=i ~ ^'-1) - T!i=i Miix) - m', the basis functions M,(x) can be canceled out and we are left with 

discrete differentiation, S, = f, - f,-] using the coboundary process as described in Section |3] The conservation laws 
thus become completely independent of the basis functions. This means that they are exactly satisfied, even on coarse 
and on arbitrary curved meshes, as already proven by the commutation property of the projection with the pullback 
(Ei. 

While the set (Ml'Lo sums up to 1 (partition of unity), the set of edge basis functions, {M, )''^[, scale such that they 
provide unit integral, 

M, := r M,(x) = 1. for / = 1,2, ...,n. (42) 

We proof this in the case that the set {A^,(x))"^q has interpolating end-conditions, i.e. T^(xo) = 2"^o Ni{xo) = Tq and 
Th{x„) = Tj"=o fi Ni(x„) = T„. Then we have by Stokes theorem (jS), 

Jxo 

Since this should hold for all {f, € R}" q. Mi - M„ - I and consequently {m, = l}"^. 

Important is to realize that the derivation of the edge basis functions is valid for any basis {Ni{x)}"^Q (polynomial 
or non-polynomial) that is linear independent and is a partition of unity. Examples are Lagrange polynomials, Bezier, 
B-splines and NURBS. If the basis is nodal, which is the case for Lagrange polynomials, then the m, have geometric 
as well as physical meaning since they represent the line integrals along the mesh edges e,. Edge functions based 
on Lagrange polynomials are used in the Mimetic framework presented in |20, 21, 22, 23J. In case of a non-nodal 
basis {Ni{x)}"^Q, such as Bezier polynomials, B-splines and NURBS, the degrees of freedom S, only have a geometric 
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interpretation, since they are discrete values associated with mesh edges e,. They are however not related to differential 
forms by means of integration. 

Compatible spaces of B-splines have already been used in HI [Tsl [13, [Si ■ They make use of the fact that the 
derivative of a B-spline of order p is a B-spline of order p - I and can be written in the following form Issll . 



(43) 



where the c, are constants that depend on p, and the regularity of the mesh. We can directly observe that using 
B-splines, the 'node' functions, Ni{x), are given by Ni{x) - Bip{x), while the 'edge' functions, Mi{x), are given 
by M,(x) = c, ■ Bip-\{x). The latter are called Curry Schoenberg B-splines, which was in fact the initial B-spline 
representation derived by Schoenberg llsill . It was later that B-splines were scaled to form a partition of unity. 

In this paper we derive compatible spaces of discrete differential forms from NURBS ll37ll . non-uniform rational 
B-splines. NURBS are rational functions of B-splines and are given as. 



Ni{x) 



Z'U^i-Bi,p(xy 



(44) 



NURBS trivially satisfy the required partition of unity, and can therefore be used as a basis to derive compatible spaces 
of discrete differential forms. Figure [TT] shows an example of cubic NURBS basis functions and their corresponding 
edge functions for different choices of the NURBS weights w,-. 

4.4. Multivariate discrete differential fonns 

A basis for differential forms in n-dimensional space is simply obtained by means of the tensor product of the 
derived univariate 'node' and 'edge' functions. A basis for 2-forms in 2D is for example obtained by applying the 
tensor product of edge functions in both directions. We can define the following finite dimensional approximation 
spaces for 0-, 1-, and 2-forms in R^, 

AO (Q) = span [mx') ® N j(x^"'j^^.^^ 

Al (Q) = span {m,(x1) ® N j{x^"'j"'.^^ x span [Ni(x') ® M j{x^"'j^^.^^ 
Al (Q) = span [Mi(x') ® M jix^"'j"-.^^ 

Examples of 2-dimensional basis functions associated with inner and outer oriented points, edges, and faces are 
displayed in Figure [12] and [13] These spaces follow the following sequence, which is exact on contractible domains. 



A°(Q) 



A°(Q) 



A\Q.) 



A2(Q) 



(45) 



ArXQ) 



We will more deeply discuss the discrete spaces of differential forms in to stay close to our discussion in section 
[1] We introduce the following notation for 3-dimensional tensor product basis functions. 



Bij,k{x) = Ni(x')Nj(x^)Ndx') 
Bf]/x) = Mi(x')Nj(x^)N,(x'), 
Bfik^^) = mx')Mj(x^)Mk{x\ 
fif//'(x) = MAx')Mj{x^)Mu{x') 



^S^(x) = Ni{x')Mj{x^)Nk{x\ 
I^l^ix) = MXx')Nj{x-)M,{x\ 



Bflfii) = Ni(x')Nj{x^)Mk{x') 
B^l^ix) = MAx')Mj{x^)Nu{x^) 
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(a) Inner oriented cell complex 





(b) Inner oriented 'node' function 



(c) Inner oriented 'node' function 














iiinmii 



(d) Inner oriented 'edge' function 



(e) Inner oriented 'edge' function 





(f) Inner oriented 'face' function 



(g) Inner oriented 'face' function 



Figure 12: 2D examples of inner oriented 'node' (b and c), 'edge' (d and e) and 'face' (f and g) basis functions. The basis functions are derived 
from NURBS of bi-degree 2 and knot vectors Ui = { 0, 0, 0, 1, 1, 1 | and U2 = { 0,0,0,0.5, 1, 1, 1 1 with weights equal to 1 (NURBS reduces to 
B-spline). 
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(a) Outer oriented cell complex 




(b) Outer oriented 'node' function (c) Outer oriented 'node' function 




(d) Outer oriented 'edge' function (e) Outer oriented 'edge' function 




(f) Outer oriented 'face' function (g) Outer oriented 'face' function 

Figure 13: 2D examples of outer oriented 'node' (b and c), 'edge' (d and e) and 'face' (f and g) basis functions. The basis functions are derived from 
NURBS of bi-degree 2 and knot vectors Ui = { 0, 0, 0, 1, 1, 1 1 U2 = { 0, 0, 0, 0.5, 1,1,1 ) with weights equal to 1 (NURBS reduces to B-spline). 
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Note that the superscripts indicate in which direction an edge function is used. We can then define the following finite 
dimensional approximation spaces for 0-, 1-, 2- and 3-forms, 



A«(Q) = span{fi,,,(x)};::;;^^^^„ 

A? (O) = span B*'',(x) " " x span B^^' (x) ' x span B'^', (x) ' 

Aj (« = span i«:s'(x)n;:;,, X span (Cf^)!:™.., « i^;!/"-^' 

A;;(£!)=spa„|B»^'(x)|';:;;;.:;,., 

These spaces are constructed such that they satisfy the following sequence, which i; 



A°(Q) 



vhich is exact on contractible domains, 

^'(") ^ ^'(") ^ ^'(") 



(46) 



^°(") ^ ^'(") ^ ^'(") ^*(") 



div 



We emphasize that 0- and 3-forms are scalar valued functions and 1- and 2-forms are vector valued functions. 
Some examples of quantities such as temperature, velocity, vorticity and density are displayed in the following exam- 
ple. 

Example 4.1 (Reconstruction in R^^). Examples of the finite dimensional reconstructions of 0-, 1-, 2- and 3-form 
fields. 



0-form: rf (x) = ^ f B,- ,;t(x) 



2- form: (X) = ^ -l^^. ^S^^) ^ Z ^ Z ^W^^ 

i,j,k i,j,k i,j,k 

3- form: >(x) = ^ B*',:^''(x) 

i,j,k 



As in univariate space, the balance laws are completely independent of the basis functions and we can perform 
differentiation discretely. We thus have a discrete matrix representation of the gradient, curl and divergence operator 
which is exact, even on curved and coarse meshes, see Section[3] These relations hold for Lagrange, Bezier, B-spline, 
NURBS or any other basis {A^/(x))"^q that forms a partition of unity. This means that these relations hold independent 
of a specific numerical method. 

Example 4.2 (Discrete gradient operator.). Consider m,'(x) = grad 7'j'(x) where m^(x) and T^^i's.) are of the form as 
in example \4.1\ Then the degrees of freedom ufjf,, where d = 1, 2, 3 represents component direction, are associated 
with edges and are given by, 

"%j,k ^ "^'J-k - Ti-\,j,k, uf^i^ = Tij^k - Tij-i,k and ufjj^ = f,j;t - Tij^k-i 

This discrete representation allows a matrix vector product u = Di o T where Dio = Eg j is an incidence matrix, see 
Examvle \3.1\ This discrete gradient is exact and does not depend on geometric transformations, as proven in ( I29D . 
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Example 4.3 (Discrete curl operator.). Consider w^(x) - curl m^(x) where w^(x) and m^(x) are of the form as in 
example \4.1\ Then the degrees of freedom oif \ are associated with faces and are given by. 





-(3) 


-(3) 


i,],k-\ 


- 2 

- u. ., 

i,j,k 




t,J,k 




i-\,],k 


-(3) 




-(2) 

- u. ., - 

ij,k 


-(2) 
i-l,],k 




i,j.k 



which allows a matrix vector product cj — D2,i u where D2,i — E^j '■^ incidence matrix, see Example \3.1\ This 
discrete curl is exact and does not depend on geometric transformations, as proven in H29\l . 

Example 4.4 (Discrete divergence operator.). Consider /)^(x) - div w^(x) where cj^ix) and //^(x) are of the form 
as in example \4.1\ Then the degrees of freedom fj^k ore associated with volumes and follow the discrete relation, 

f _ ,-,(1) ,-,(1) , ,-,(2) -(2) - (3) - (3) 

•I'd^'^ uj,k i-i~J~k uj,k i,J-^,k i,],k ij,k-l 

which allows a matrix vector product f — D3 2 cj where D2 2 — is an incidence matrix, see Example \3.1\ This 
discrete divergence is exact and does not depend on geometric transformations, as proven in \29h 

To summarize, with the appropriate basis functions in which to represent continuous differentiation, yS^(x) = 
d a^j"'(x) {k — 1,2, 3), the gradient, curl and divergence reduce to a relation between the expansion coefficients. 



5. Mixed formulation for Stokes flow in terms of differential forms 

In this section we will use the finite dimensional spaces of differential forms, developed in the previous section, 
in a mixed continuous Galerkin setting. We will apply the method to several 2D test problems concerning Stokes 
flow (see section [1). Which mixed formulation we choose to attack the problem, will solely dependent on physical 
considerations. Following the reasoning in section [1] we choose all variables as outer oriented differential forms, 
such that we can apply the topological divergence for mass conservation, see Figure [3] This ultimately results in a 
point-wise divergence free velocity field in the discrete setting. Furthermore, choosing all variables as outer oriented 
differential forms, has the advantage that normal velocity is prescribed in a strong sense, while tangential velocity is 
prescribed in a weak sense. 

Consider an «-dimensional domain Q., filled with an incompressible fluid with constant viscosity v. Under these 
assumptions. Stokes flow can be described in terms of differential forms as, 

d*u"-^ = cj"-', d oj"-^ + d* p" = f"-\ d u"-^ = (47) 

Here velocity is an outer oriented (n - l)-form, m""', pressure is an outer oriented n-form, p", the vorticity is an outer 
oriented (« - 2)-form and the right-hand-side is an outer oriented (n - l)-form, /' 



n-l 



We follow the same reasoning as in 112 lH and derive the mixed formulation of the Stokes problem in 2 steps, 

1. Multiply the equations in ( |47] ) by test functions a""~, y6""' and y" to obtain inner products terms. 

2. Use integration by parts to replace the coderivative d* by the exterior derivative d and an additional boundary 
integral using ( fT6b . 

The Stokes problem can now be posed as foflows: find \cj"-^ e A""^, m""' e A""', e A"| for all \a"-^ e A""^ , 

pn-\ g ^,,-1^ y, g g^^jj jjj^j 

A (48) 



{i3"-\ d u"-\ + (d /3"-\ p'% - r A = [r\ r'\ 

{y",du"-\^0. (50) 



an 

(49) 
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(a) grid 1 (b) grid 2 

Figure 14: The Cartesian and curved mesh used in the convergence study. 



This is the vorticity-velocity -pressure (VVP) formulation also used in lll7ll2lLl22L 13811 which is well posed 11221 13911 . 
In the discrete problem we merely replace the infinite dimensional spaces of differential forms by (wjp^ e K"^^ , 
p"^ e h'l\. Because, Ah(Q.) c A(Q), the discrete problem is automatically well posed as well. The 



„«-i 



problem is closed by imposing normal and tangential velocity components at the boundary. See 11211 12211 for an 
overview of all admissible boundary conditions. The resulting system is symmetric and given by. 



• -M„_2 
M„_iD„_i_„_2 



ri-lM-2 





M„D„,„_i 





D[„_,M„ - B2(p) 














u 








p. 




I j 



(51) 



The discrete curl D„_i,„_2 = E^_2n-i ^'^'^ divergence D„ „_i - j ^ are incidence matrices, see Section[3] depending 
solely on mesh topology and are exact irrespective of the geometry mapping <S> and coarseness of the mesh. The inner 
product matrices M„_2, M„_i and M„ are calculated by pulling the finite dimensional spaces of diff'erential forms 
back to the reference domain Q' (section|2]i and performing Gauss numerical integration there. Bi(m""') and B2(p) 
represent the natural boundary terms, tangential velocity and tangential pressure, that arise due to integration by parts 
using ( fTSI l. Since the tangential velocity is weakly enforced, it appears in the right hand side, while the pressure term 
B2 is not enforced and thus remains in the left hand side. 



6. Numerical results 

We will test the method on several common numerical test cases: 1) We compare with an analytical solution (same 
test problem as in 12111 ') and show results for /i-refinement on a Cartesian and skewed, non-linear mesh; 2) Taylor 
Couette flow and compare with the analytical solution; and 3) lid driven cavity flow and compare with benchmark 
results of 

6.1. Manufactured solution 

Consider as computational domain the unit square, Q = [0, 1]^, filled with an incompressible fluid of viscosity, 
V = 1 . As an analytical solution we choose, 

tdPix, y) - - 4-n sin (2nx) sin (Iny). 

u^(x,y) - - sin (2nx) cos (2;ry) ■ dy - cos (Ittx) sin (2ny) ■ dx 
p'(x, y) — sin (nx) sin (ny) ■ dx dy. 

f^ix,y) = (-Stt^ sin (Inx) cos (Iny) + n cos (ttx) sin (jiy)^ ■ dy+ 
(-Stt^ cos (Inx) sin (Iny) - n sin {nx) cos (ttj)) ■ dx. 

We perform calculations on a Cartesian mesh as well as on a bicubic degree curved mesh, see Figure [141 
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(a) Vorticity - grid 1 



(b) Velocity - grid 1 



(c) Pressure - grid 1 
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(d) Vorticity - grid 2 
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(e) Velocity - grid 2 
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(f) Pressure - grid 2 



Figure 15: H-convergence for vorticity, velocity and pressure on grid 1 (a,b,c) and 2 (d,e,f). 



The divergence of m^' is point wise zero at all stages of refinement. Figure [15] shows the convergence behavior 
under /j-refinement of the vorticity, velocity and pressure for polynomial orders ranging from p = 2 to p = 5 for the 
vorticity and p -\\a p - \ for the velocity and pressure. The maximum mesh-size hmax is calculated as the maximum 
diagonal length of all elements divided by the square root of 2. We can observe optimal order of convergence for all 
variables on both grids. Only the pressure on grid 2 of polynomial order p — \ shows some odd behavior. This is 
probably due to the low order of approximation of the variables in conjunction with a higher order {p - 4) non linear 
mapping of the geometry (super-parametric). 

6.2. Taylor Couette flow 

We consider a somewhat more realistic problem with an incompressible fluid with viscosity, v = 1, in a domain 
bounded by two contra rotating cylinders. The outer cylinder with radius r,,,,, = 2 is fixed, while the inner cylinder 
with radius r,„ - 1 rotates in counter clockwise direction with an angular velocity equal to 1 . 

The flow field in Cartesian coordinates (x,y) as a function of radius r is then described by the following velocity 
field, 

u\r) = i--r+ —]{-sin0dy + cos0dx) (52) 
\ 3 3r/ 

and is depicted in Figure [T6k. 

The geometry is exactly represented by 4 C" NURBS patches of degree 2 by 1 . Note that this test-case is more 
interesting since not all weights are equal to one. The weights associated with the four inner and four outer corners of 
the domain have weights equal to 0.5 ■ V2. Figure [T6b shows optimal convergence of the velocity under /i-refinement 
for polynomial orders of p = 1 to p = 4 in the radial direction. The pressure is constant throughout the domain and is 
up to machine precision at all stages of refinement. Here h,„ax represents the mesh-size in radial direction. 
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(a) velocity field and NURBS geometiy (b) ft-convergence velocity 

Figure 16: Flow field and NURBS geometiy representation (a) using 4 patches; and /i-convergence for velocity (b). 




(a) streamfunction (b) vorticity (c) pressure 



Figure 17: Results lid driven cavity flow on a 60x60 uniform bi-cubic NURBS mesh (weights ai'e all equal to one) of highest regularity. 
6.3. Lid driven cavity flow 

The lid driven cavity flow is one of the classical benchmark cases to asses numerical methods and to verify Navier 
Stokes codes. Consider the unit square domain Q. - [0, 1]^ filled with an incompressible fluid of viscosity v - I. On 
the top of the domain we apply a tangential velocity u{x, 1) = 1 while on the other sides of the domain the tangential 
velocity is set to zero. The result is a clockwise rotating flow with small counter rotating eddies in the two lower 
corners. Because of the discontinuity of the velocity in the two upper corners, both the vorticity and pressure are 
infinite at these places. These singularities make the lid driven cavity flow a challenging test case. 

Figure [T7] shows results for the stream function, vorticity and pressure on a bi-cubic uniform grid of maximum 
regularity and 60x60 degrees of freedom. We note that the divergence of velocity is point-wise zero in whole Q. 
Furthermore no special treatment has been given to the corner singularities. 

In Figure [18] the horizontal component of velocity has been plotted along the vertical centerline (0.5,}') and the 
vertical component of velocity has been plotted along the horizontal centerline {y, 0.5). The results confirm very well 
with the benchmark results of ll40ll even though the mesh is very coarse (9x9 uniform grid of polynomial order 1, 3 
and 5). The most striking result is, however, that of the low order approximation. Note that although there is a clear 
point-wise difference, the integral values seem to match very well. This is a direct consequence of the conservation 
properties that we have build into the basis. 
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(a) (b) (c) 



9x9, p=l 9x9, p= 3 9x9, p= 5 




(d) (e) (f) 

Figure 18: Comparison of numerical approximation (blue) with benchmark results of (icil (red). Horizontal (a,b,c) and vertical (d,e,f) velocity 
profile at centerlines of cavity for a 9x9 uniform grid of order 1, 3 and 5. 



7. Conclusion 

We have developed arbitrary order interpolants, from any basis that is a partition of unity, that satisfies a discrete 
Stokes theorem. The resuhing gradient, curl and divergence conforming spaces have the property that the conservation 
laws become completely independent of the basis functions. As an example, we have derived conforming spaces form 
NURBS, and thereby generalized the discrete spaces of differential forms introduced in ITf). We have applied these 
new spaces in a mixed Galerkin setting which amongst others resulted in an exactly divergence free discretization of 
Stokes flow. 
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